function [ results_fnames ] = osl_run_first_level_continuous( oat )

% [ results_fnames ] = osl_run_first_level_continuous( oat )
%
% takes in an OAT, which needs to be setup by calling oat=osl_setup_oat(S), struct
% and runs first level continuous time GLM or seed based connectivity
%
% This function should normally be called using osl_run_oat(oat);
%
% MWW 2012

global OSLDIR;

first_level=oat.first_level;
source_recon_name='source_recon';

% check if user has asked not to do GLM
if(~isfield(first_level,'doGLM')),
    first_level.doGLM = 1;
end

if(first_level.time_average)
    error('Time averaging over whole window is not appropriate on continuous data');   
end;

% set the downsample factor.  If doing a tf, read tf_downsample_factor.  If
% doing time domain, read time_downsample_factor
switch first_level.tf_method
    case 'none' % use time_downsample_factor
        error('Continuous data time-wise GLM only handles time-freq analysis');
        
        if ~isempty(first_level.time_downsample_factor)
            ds_factor = 1./first_level.time_downsample_factor;
        else
            ds_factor = [];
        end
    otherwise  % use tf_downsample_factor
        if ~isempty(first_level.tf_downsample_factor)
            ds_factor = 1./first_level.tf_downsample_factor;
        else
            ds_factor = [];
        end
end % switch first_level.tf_method

if(~isfield(first_level,'connectivity_seed_mni_coord')),
    first_level.connectivity_seed_mni_coord=[];
end;

Ntrials   = 1;    
   
clear stats;

for subi_todo=1:length(first_level.sessions_to_do),   
        
    sub=first_level.sessions_to_do(subi_todo);
    
    %sub/length(first_level.submatfiles),
    
    %% load in the beamformer result for this session
    source_recon_results_fname = oat.source_recon.results_fnames{sub};
    source_recon_results=osl_load_oat_results(oat,source_recon_results_fname);
    
    % get D and update path in case OAT dir has been moved
    D=source_recon_results.BF.data.D;
    if(~strcmp([oat.source_recon.dirname],D.path)),        
        D=path(D,[oat.source_recon.dirname]);
        D.save;
        D=spm_eeg_load([D.path '/' D.fname]);
        osl_save_oat_results(oat,source_recon_results);    
    end;
    
    %% setup things that are common to all sessions
    if(subi_todo==1),
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% set up mask indices from mask/coords
        if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
            
            disp('Working in sensor space');
  
            %%%%%%%%%%%%%%%
            %% setup channels
            modality='MEG';
            chanindmeg = strmatch(modality, D.chantype); % need to even include bad channels in output so that it is the same for all sessions
       
            first_level_results.mask_indices_in_source_recon=1:length(chanindmeg);
    
        else,
           
            % setup std space brain
            first_level_results.gridstep=source_recon_results.gridstep;
            
            S=[];
            S.lower_level_mask_fname=[oat.source_recon.dirname '/' source_recon_name '_mask'];
            S.current_level_mask_fname=[oat.source_recon.dirname '/' oat.first_level.name '_mask'];
            S.current_level=first_level;
            S.lower_level_mni_coord=source_recon_results.mni_coord;
            S.lower_level_gridstep=source_recon_results.gridstep;
            
            [first_level_results.mask_indices_in_source_recon, first_level_results.mni_coord]=setup_mask_indices(S);
            clear S;
        end;
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% set up potential time frequency bases in tf_settings, note that any time or
        % frequency averaging will be done AFTER the time-frequency decomposition
        
        % set the frequency range
        freq_range = first_level.tf_freq_range; % range used for stat testing
        if isempty(freq_range)
            % if not set as a first_level option, default to the range used
            % for source_recon
            freq_range=source_recon_results.source_recon.freq_range;
            
            if isempty(freq_range)
                error('No frequency range has been specified. Set oat.source_recon.freq_range or oat.first_level.tf_freq_range');
            end % isempty(freq_range)
        end % if isempty(freq_range)
        
        source_recon_time_indices=find(source_recon_results.samples2use);
        source_recon_times = D.time(source_recon_time_indices);
        
        if isempty(first_level.time_range)
            first_level.time_range=source_recon_results.woi;
        end
        
        tf_settings = [];
        tf_settings.tf_method               = first_level.tf_method;
        tf_settings.tf_logtransform         = first_level.tf_logtransform;
        tf_settings.tf_hilbert_freq_res     = first_level.tf_hilbert_freq_res;
        tf_settings.tf_freq_range           = freq_range; % set above
        tf_settings.tf_num_freqs            = first_level.tf_num_freqs;
        tf_settings.raw_times               = source_recon_times;
        tf_settings.ds_factor               = 1; % do not downsample inside transform - instead we will do this later after the moving average low pass filtering
        tf_settings.tf_morlet_factor        = first_level.tf_morlet_factor;
        tf_settings.tf_multitaper_taper                = first_level.tf_multitaper_taper;
        tf_settings.tf_multitaper_twin                 = first_level.tf_multitaper_twin;
        tf_settings.tf_multitaper_freqsmooth           = first_level.tf_multitaper_freqsmooth;
        
        out = osl_tf_transform( tf_settings ); % passed without a data matrix, osl_tf_transform just gets container sizes
        
        if isfield(out, 'tf_morlet_basis')
            first_level.tf_morlet_basis = out.tf_morlet_basis;
        end
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% setup time and frequency windows/indices for using later
          
        % setup result time dimension        
        first_level_time_indices  = ( out.tf_times >= first_level.time_range(1) ) & ( out.tf_times <= first_level.time_range(2) );
        first_level_results.times=1;       
        if first_level.time_average
            disp('Note time averaging irrelevant and so is ignored for continuous data time-wise GLM');
        end
        
        % setup result freq dimension
        first_level_results.frequencies = out.tf_freqs;
        if(first_level.freq_average),
            first_level_results.frequencies=mean(first_level_results.frequencies);
        end;
       
        % setup result spatial dimension
        Nvoxels_out=length(first_level_results.mask_indices_in_source_recon);
        if(first_level.space_average),
            disp('Doing spatial average over all mni coords');
            Nvoxels_out=1;
            if(isfield(first_level_results,'mni_coord')),
                first_level_results.mni_coord=mean(first_level_results.mni_coord,1);
            end;
        end;
        
        clear out

        %%%%%%%%%%%%%%%%
        %% setup seed for seed-based connectivity analysis       
        if(~isempty(first_level.connectivity_seed_mni_coord)),
            
            if(first_level.tf_num_freqs>1)
                error('Continuous data seed-based connectivity only handles single freq bands at a time');
            end;

            mni_coords=first_level_results.mni_coord;
            
            disp(['Using for seed: MNI coordinate ' num2str(first_level.connectivity_seed_mni_coord)]);
            
            % find the nearest index of the beamformed voxels to the specified mni_coord
            dists=(sqrt(sum((mni_coords-repmat(first_level.connectivity_seed_mni_coord,size(mni_coords,1),1)).^2,2)));
            
            [dist,first_level_results.seed_index]=min(dists);
            
            disp(['Nearest dipole for seed at MNI coordinate ' num2str(mni_coords(first_level_results.seed_index,:))]);                                              
                       
            first_level.contrast=[];
            first_level.contrast{1}=[1]';

        end;        
        
        %%%%%%%%%%%%%%%%
        %% results containers
        ntpts=1;
        nfreqs=length(first_level_results.frequencies);
        
        if(first_level.doGLM)
            first_level_results.stdcope=zeros([Nvoxels_out,ntpts,length(first_level.contrast),nfreqs]);% num_voxels x  num_timepoints x num_contrasts x num_freqs x Nvoxels_multidipole
            first_level_results.cope=zeros([Nvoxels_out,ntpts,length(first_level.contrast),nfreqs]);% num_voxels x num_timepoints x num_contrasts x num_freqs x Nvoxels_multidipole
        end;

        first_level_results.pseudo_zstat_var=zeros([Nvoxels_out,1]);
        first_level_results.pseudo_zstat_var_nai=zeros([Nvoxels_out,1]);
        
    end; % end of if(sub==1)
    
    disp(['Reconstruct time courses and computing stats for dataset ' source_recon_results.fname]);    

    %%%%%%%%%%%%%%%
    %% setup channels
    modality='MEG';
    chanindmeg = strmatch(modality, D.chantype);
    chanind = setdiff(chanindmeg, D.badchannels);
    if isempty(chanind)
        error(['No good ' modality ' channels were found.']);
    end
    if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
        % for sensor space - also include bad sensors sensors 
        first_level_results.mask_indices_in_source_recon=1:length(chanindmeg);
        chanind = chanindmeg;
    end;    
    
    %%%%%%%%%%%%%%%
    %% load in sensor data
    first_level_results.D_sensor_data=D;
    first_level_results.chanind=chanind;
    
    sensor_data=D(chanind, source_recon_time_indices, 1);
    
    if strcmp(source_recon_results.recon_method,'none'),
        NK=1;
    else
        NK=numel(source_recon_results.BF.inverse.W.MEG); 
    end;
    
    if NK>1
        classchanind=find(strcmp(D.chanlabels,'Class'));
        if(isempty(classchanind)),
            error(['No ''Class'' chanlabel in: ' D.fname]);
        else
            for kk=1:NK,        
                class_samples_inds{kk} = (D(classchanind,source_recon_time_indices, 1)==kk);
            end;   
        end;
    else,        
        class_samples_inds{1}= ones(1,length(source_recon_time_indices),1);
    end;
    
    %%%%%%%%%%%%%%%%
    %% setup the GLM design matrix
    % design matrix should be num_PEs x num_timepoints    
    if(~isempty(first_level.connectivity_seed_mni_coord)),

        % set the design matrix regressor to be the seed voxel power time course
        
        % recon seed time course
        S2=[];
        S2.source_recon_results=source_recon_results;
        S2.chanind=chanind;
        S2.class_samples_inds=class_samples_inds;
        S2.voxind=first_level_results.seed_index;
        S2.sensor_data=sensor_data;
        S2.use_class_probs=0;
        
        [ dat_seed wnorms_seed wnorm_seed wnorms_nai_seed wnorm_nai_seed] = get_voxel_recon_timecourse( S2 );
       
        %dat_seed=dat_seed'./sqrt(wnorms_nai_seed');
        %wnorms_seed=wnorms_seed'./wnorms_nai_seed';
        dat_seed=dat_seed';
        wnorms_seed=wnorms_nai_seed';
        
        %% Transform the data - either downsampling in time domain, or a TF transform
            
        tf_settings = [];
        tf_settings.tf_method               = first_level.tf_method;
        tf_settings.tf_logtransform         = first_level.tf_logtransform;
        tf_settings.tf_hilbert_freq_res     = first_level.tf_hilbert_freq_res;
        tf_settings.tf_freq_range           = freq_range; % set above
        tf_settings.tf_num_freqs            = first_level.tf_num_freqs;
        tf_settings.raw_times               = source_recon_times;
        tf_settings.ds_factor               = 1; % do not downsample inside transform - instead we will do this later after the moving average low pass filtering
        tf_settings.tf_morlet_factor        = first_level.tf_morlet_factor;
        tf_settings.tf_multitaper_taper                = first_level.tf_multitaper_taper;
        tf_settings.tf_multitaper_twin                 = first_level.tf_multitaper_twin;
        tf_settings.tf_multitaper_freqsmooth           = first_level.tf_multitaper_freqsmooth;

        if isfield(first_level,'tf_morlet_basis')
            tf_settings.tf_morlet_basis = first_level.tf_morlet_basis;
        end

        out = osl_tf_transform( tf_settings , dat_seed' ); % do the transformation
       
        dat_seed_tf = permute(out.dattf,[2 3 1]);
            
        wnorms_seed_tf=wnorms_seed;
        
        first_level.design_matrix=dat_seed_tf';

        % add a regressor to model the mean and variance normalise
        % the seed power time course
        %processed_design_matrix=[(first_level.design_matrix); ones(size(first_level.design_matrix,2),1)'];  
        processed_design_matrix=[(first_level.design_matrix)];
              
        %first_level.contrast=[];
        %first_level.contrast{1}=[1 0]';
        
    else,
        
        if(~isfield(first_level,'design_matrix')),
                
            if(first_level.doGLM)
                error('Seed voxel or design matrix needs to be specified');
            else
                % setup dummy design matrix
                first_level.design_matrix=ones(1,length(source_recon_time_indices));
            end;
        end;
                        
        % apply time windowing to design matrix
        wnorms_seed_tf=ones(size(first_level.design_matrix,2),1);
        wnorms_seed_tf=wnorms_seed_tf(source_recon_time_indices);
        
        processed_design_matrix=first_level.design_matrix(:,source_recon_time_indices);
    
    end;     
       
    if(size(processed_design_matrix)~=length(first_level_time_indices)),
        error(['Design matrix dimensions [npes x ntps] = ' num2str(size(first_level.design_matrix)) ' do not match the data: ntpts=' num2str(length(first_level_time_indices))]); 
    end;
    
    for pp=1:size(processed_design_matrix,1),
        processed_design_matrix2(pp,:)=processed_design_matrix(pp,first_level_time_indices);
    end;   
    processed_design_matrix=processed_design_matrix2;
    
    wnorms_seed_tf=wnorms_seed_tf(first_level_time_indices);    
            
    % apply moving average to design matrix
    if ~isempty(first_level.time_moving_av_win_size)
        fsample=1/(source_recon_times(2)-source_recon_times(1));        
        windowSize = first_level.time_moving_av_win_size*fsample;
                        
        wnorms_seed_tf=moving(wnorms_seed_tf,round(windowSize));
        
        for pp=1:size(processed_design_matrix,1),
            processed_design_matrix(pp,:)=moving(processed_design_matrix(pp,:),round(windowSize));
        end;        
    end;

    % temporally downsample                                     
    if ~isempty(first_level.tf_downsample_factor); 
        processed_design_matrix = osl_spm_resample(processed_design_matrix,ds_factor);
        wnorms_seed_tf = osl_spm_resample(wnorms_seed_tf',ds_factor)';

    end;

    if(oat.first_level.do_glm_demean)
        %processed_design_matrix=demean(processed_design_matrix);
        %processed_design_matrix=demean(processed_design_matrix')';
    end;
    
    % precompute some things
    x=processed_design_matrix';
    pinvxtx=pinv(x'*x);
    pinvx=pinv(x);

    contrast_list=first_level.contrast;
    for c=1:length(contrast_list),
        contrast_list{c}=contrast_list{c}(:);
    end;

    %%%%%%%%%%%%%%%%
    %% do everything from now on one voxel at a time (for the sake of
    %% RAM) in the first level mask space                
    if(first_level.doGLM)
        disp('First level COPEs outputted will have dimension Nvoxels x Ntpts x Ncontrasts x Nfreqs:');
        tmp=[size(first_level_results.cope) 1 1];
        disp([num2str(tmp(1:4))]);
    end;
         
    %%%%%%%%%%%%%%%%
    %% now do space averaging 
    if(first_level.space_average),
        error('space_average Not implemented');
        
        dat=zeros(Ntrials,length(source_recon_time_indices));
        wnorms=zeros(Ntrials,length(source_recon_time_indices));            
        num=0;
        
        for indind=1:length(first_level_results.mask_indices_in_source_recon), % indexes brain space
            
            weights_usable=0;
            
            S2=[];
            S2.source_recon_results=source_recon_results;
            S2.chanind=chanind;
            S2.class_samples_inds=class_samples_inds;
            S2.voxind=first_level_results.mask_indices_in_source_recon(indind);
            S2.sensor_data=sensor_data;
            S2.use_class_probs=0;
            
            [ datnew wnormsnew wnorm wnorm_nais wnorm_nai ] = get_voxel_recon_timecourse( S2 );
            
            if sum(isnan(squash(dat)))==0,
                num=num+1;
                dat=dat+datnew;
                wnorms=wnorms+wnormsnew;
            end;
            
        end;
        
        dat=dat/num;
        wnorms=wnorms/num;
    end;
      
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% now loop over voxels
    ft_progress('init', 'etf');
    
    % debug:
    %[indfind, mni_coord]=nearest_vec(source_recon_results.mni_coord,[-32 -26 50]);
                            
    for indind=1:Nvoxels_out, % indexes brain space      
        
        ft_progress(indind/Nvoxels_out);                                   

        if(~first_level.space_average),
            S2=[];
            S2.source_recon_results=source_recon_results;
            S2.chanind=chanind;
            S2.class_samples_inds=class_samples_inds;
            S2.voxind=first_level_results.mask_indices_in_source_recon(indind);
            S2.sensor_data=sensor_data;
            S2.use_class_probs=0;
            
            [ dat wnorms wnorm wnorms_nai wnorm_nai] = get_voxel_recon_timecourse( S2 );
       
            %dat=dat'./sqrt(wnorms_nai');
            %wnorms=wnorms'./wnorms_nai';
            
            dat=dat';
            wnorms=wnorms_nai';
        
            first_level_results.pseudo_zstat_var(indind)=var(dat./sqrt(wnorms));
            %first_level_results.pseudo_zstat_var_nai(indind)=var(dat./sqrt(wnorms_nai));
           
        end;
         
        %dat=dat';
        %wnorms=wnorms';
        
        if sum(isnan(squash(dat)))~=0,

            if(first_level.doGLM)
                first_level_results.cope(indind,:,:,:)=nan;
                first_level_results.stdcope(indind,:,:,:)=nan;         
            end;
            warning('NANs is recon data');
        else
           
            if(~isempty(first_level.connectivity_seed_mni_coord)),
                if(first_level.connectivity_seed_regress_zerolag),

                    % orthoganalise the target voxel WRT the seed voxel (doing
                    % it this way round means that the TF transform and
                    % windowing does not need to be redone for the seed voxel)
                    if(NK>1),
                        hmm_states=source_recon_results.block(1).q_star;
                        sts=unique(hmm_states);
                        for i=1:length(sts),
                            tinds=(hmm_states==sts(i));
                            dat1=dat(tinds);
                            dat2=dat_seed(tinds);

                            dat(tinds)=dat1-dat2*(dat2\dat1);
                        end;     
                    else,
                        dat1=dat;
                        dat2=dat_seed;

                        dat=dat1-dat2*(dat2\dat1);
                    end;        
                    
                end;
            end;
            
            %% Transform the data - either downsampling in time domain, or a TF transform

            tf_settings = [];
            tf_settings.tf_method               = first_level.tf_method;
            tf_settings.tf_logtransform         = first_level.tf_logtransform;
            tf_settings.tf_hilbert_freq_res     = first_level.tf_hilbert_freq_res;
            tf_settings.tf_freq_range           = freq_range; % set above
            tf_settings.tf_num_freqs            = first_level.tf_num_freqs;
            tf_settings.raw_times               = source_recon_times;
            tf_settings.ds_factor               = 1; % do not downsample inside transform - instead we will do this later after the moving average low pass filtering
            tf_settings.tf_morlet_factor        = first_level.tf_morlet_factor;
            tf_settings.tf_multitaper_taper                = first_level.tf_multitaper_taper;
            tf_settings.tf_multitaper_twin                 = first_level.tf_multitaper_twin;
            tf_settings.tf_multitaper_freqsmooth           = first_level.tf_multitaper_freqsmooth;

            if isfield(first_level,'tf_morlet_basis')
                tf_settings.tf_morlet_basis = first_level.tf_morlet_basis;
            end

            out = osl_tf_transform( tf_settings , dat' ); % do the transformation

            dat_tf = permute(out.dattf,[2 3 1]);
            wnorms_tf=wnorms;

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%            
            %% now do windowing over time 
            dat=dat_tf(first_level_time_indices,:);                
            wnorms_tf=wnorms_tf(first_level_time_indices,:);    
            clear wnorms;
           
            if(first_level.freq_average)
                dat=mean(dat,2);                
            end;

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
            %% moving average within
            %% the window specified
            if ~isempty(first_level.time_moving_av_win_size)

                fsample=1/(source_recon_times(2)-source_recon_times(1));        
                windowSize = first_level.time_moving_av_win_size*fsample;

                dat=moving(dat,round(windowSize));
            
                wnorms_tf=moving(wnorms_tf,round(windowSize));
                
            end

            % temporally downsample                                     
            if ~isempty(first_level.tf_downsample_factor);
                dat= osl_spm_resample(dat',ds_factor)';  
                wnorms_tf= osl_spm_resample(wnorms_tf',ds_factor)';
               
            end;
            glm_input_times=linspace(first_level.time_range(1),first_level.time_range(2),length(dat));           
                                                                         
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
            %% now do the GLM stats
            if first_level.doGLM
                
                %if(oat.first_level.do_glm_demean)
                %    dat=demean(dat);
                %end;
                
                if(size(dat,1)>1)                    
                    for f = 1 : size(dat,2),
                        t=1;
                        
                        dattf=dat(:,t,f);
                        x=processed_design_matrix';
                        
                        if(0),%indind==indfind)
                            keyboard;
                            x2=(x);
                            dattf2=(dattf);
                            
                            figure;plot((x2));ho;plot((dattf2),'r');                            
                            [cope,vc]=ols((dattf2), (x2), [1]);
                            cope/sqrt(vc)
                            
                            %%
                            figure;plot(wnorms_tf(:,t,f));ho;plot(wnorms_seed_tf(:,t,f),'r');
                            dattf2=dattf./sqrt(10^5*wnorms_seed_tf(:,t,f));
                            x2=x./sqrt(10^5*wnorms_seed_tf(:,t,f));
                            figure;plot(demean(x2));ho;plot(demean(dattf2),'r'); 
                            
                            [cope,vc]=ols(demean(dattf2), demean(x2), [1]);
                            cope/sqrt(vc)
                        end;
                        
                        if(NK>1),
                            
                            % apply normalisation
                            if(~isempty(first_level.connectivity_seed_mni_coord)), 
                                dattf=dattf;%./sqrt(wnorms_tf(:,t,f)+wnorms_seed_tf(:,t,f));
                            else
                                dattf=dattf;%./sqrt(wnorms_tf(:,t,f));
                            end;                           
                                                        
                            for pp=1:size(x,2),                                
                                
                                if(~isempty(first_level.connectivity_seed_mni_coord)),                
                                    x(:,pp)=x(:,pp);%./sqrt(wnorms_tf(:,t,f)+wnorms_seed_tf(:,t,f));
                                else
                                    x(:,pp)=x(:,pp);%./sqrt(wnorms_tf(:,t,f));
                                end;                                
                                
                            end;
                            
                        end;
                        
                        if(oat.first_level.do_glm_demean)
                            dattf=demean(dattf);                                                        
                        end;
                        
                        if(oat.first_level.do_glm_demean)
                            for pp=1:size(x,2),
                                x(:,pp)=demean(x(:,pp));
                            end;
                        end;
                        
                        pinvxtx=pinv(x'*x);
                        pinvx=pinv(x);
                        
                        if(first_level.use_robust_glm),
                            [b,stats2]=robustfit(x,dattf,'bisquare',4.685,'off'); % fits GLM robustly, Y=XB+E
                            for c=1:length(contrast_list),

                                if(strcmp(Sin.cope_type,'coape')),
                                    first_level_results.cope(indind,t,c,f)=contrast_list{c}'*abs(b);
                                elseif(strcmp(Sin.cope_type,'cope')),
                                    first_level_results.cope(indind,t,c,f)=contrast_list{c}'*b;
                                elseif(strcmp(Sin.cope_type,'acope')),
                                    first_level_results.cope(indind,t,c,f)=abs(contrast_list{c}'*b);
                                else
                                    error('Invalid cope_type');
                                end;
                                first_level_results.stdcope(indind,t,c,f)=sqrt(contrast_list{c}'*stats2.covb*contrast_list{c});
                            end;
                            pe(t,:)=b; % num_timepoints x  num_pes
                        else
                            [copeout, varcopeout, coapeout, dof,peout]=glm_fast_for_meg(dattf,x,pinvxtx,pinvx,contrast_list,0);
                            if(strcmp(oat.first_level.cope_type,'coape')),
                                first_level_results.cope(indind,t,:,f)=coapeout;
                            elseif(strcmp(oat.first_level.cope_type,'cope')),
                                first_level_results.cope(indind,t,:,f)=copeout;
                            elseif(strcmp(oat.first_level.cope_type,'acope')),
                                first_level_results.cope(indind,t,:,f)=abs(copeout);
                            else
                                error('Invalid cope_type');
                            end;

                            first_level_results.stdcope(indind,t,:,f)=sqrt(varcopeout);
                            pe(t,:)=peout; % num_timepoints x  num_pes
                        end;

                        %% do bc on the copes
                        if(sum(first_level.bc)>0),
                            error('No first_level.bc available for continuous data');                            
                        end;
                    end;


                else
                    first_level_results.cope(indind,:,1,:)=dat(:,:);
                    first_level_results.stdcope(indind,:,1,:)=1;
                end;
            
            else,
                if(~isfield(first_level_results,'glm_input_data')),
                    first_level_results.glm_input_times=glm_input_times;
                    first_level_results.glm_input_frequencies=out.tf_freqs;
                    first_level_results.glm_input_data=zeros([Nvoxels_out, size(dat)]);% num_voxels x  num_timepoints x num_contrasts x num_freqs                    
                end;                    
                    
                first_level_results.glm_input_data(indind,:,:,:)=dat;
            end;                    
        end;
    end;
    
    if (first_level.doGLM && ~isempty(first_level.connectivity_seed_mni_coord)),
        first_level_results.cope(first_level_results.seed_index,:,:,:)=1;
        first_level_results.stdcope(first_level_results.seed_index,:,:,:)=1e16;
    end;
    
    ft_progress('close');
    
    %% save results
    if first_level.doGLM
        first_level_results.x=x;
        first_level_results.contrasts=first_level.contrast;
    end;
    
    first_level_results.source_recon_results_fname=source_recon_results_fname;
    first_level_results.first_level=first_level;
    first_level_results.source_recon=source_recon_results.source_recon;
    first_level_results.level=1;
    %first_level_results.source_recon_results_fname=source_recon_results.fname;
    first_level_results.name=first_level.name;
    first_level_results.recon_method=source_recon_results.recon_method;
    first_level_results.session_name=source_recon_results.session_name;
          
    first_level_results.fname=[ first_level_results.session_name '_' first_level_results.name ];
    disp(['Saving statistics in file ' oat.source_recon.dirname '/' first_level_results.fname]);
    osl_save_oat_results( oat, first_level_results );
    results_fnames{sub}=first_level_results.fname;

    if (strcmp(source_recon_results.recon_method,'none') | first_level.space_average),
    else
        disp(['To create niftii files from this use a call to osl_save_nii_stats']);
    end;
    
end;

